! -------------------------------------------------
!> 八节点六面体单元,完全积分,B-Bar技术修正体积自锁
!> 节点自由度ux,uy,uz;
!> 单元节点排序和abaqus的C3d8一致;
!> 应力分量排序为 {11,22,33,12,23,13} ,和abaqus/explicit一致
!> 单元定义需: 节点坐标,模量E,泊松比nu,密度rho,热膨胀系数alpha
!> 支持的单元载荷: 体积力(N/m^3);压力载荷(N/m^2);表面力载荷(N/m^2)
!> 支持集中节点力加载:F1,F2,F3
!> 节点力向量,节点位移向量使用1维数组
!> 单元的等参坐标使用r,s,t表示
! 单元面的编号,按照右手定则会得到向内的法向
!               Face 1: 1 – 2 – 3 – 4 face
!               Face 2: 5 – 8 – 7 – 6 face
!               Face 3: 1 – 5 – 6 – 2 face
!               Face 4: 2 – 6 – 7 – 3 face
!               Face 5: 3 – 7 – 8 – 4 face
!               Face 6: 4 – 8 – 5 – 1 face
! -------------------------------------------------

!> 该模块包含了C3D8I的单元实现
module solid_mod
    use iso_fortran_env,only:real64,int64
    implicit none
    !> 单元的节点数,一个单元8个节点
    integer,parameter :: nn_el=8
    !> 一个节点的自由度数目: ux,uy,uz
    integer,parameter :: nndofs=3
    !>一个单元的自由度数目
    integer,parameter :: eldofs=24
    !> 单元的vtk代码
    integer,parameter :: c3d8i_vtk_code = 12 
    !> 单元节点i的等参坐标r数组
    real(real64),private,parameter :: r_coords(nn_el) = [-1.0_real64,  1.0_real64,  1.0_real64,&
                                    -1.0_real64, -1.0_real64,  1.0_real64, 1.0_real64, -1.0_real64]
    !> 单元节点i的等参坐标s数组
    real(real64),private,parameter :: s_coords(nn_el) = [-1.0_real64, -1.0_real64,  1.0_real64,&
                                    1.0_real64, -1.0_real64, -1.0_real64, 1.0_real64,  1.0_real64]
    !> 单元节点i的等参坐标t数组
    real(real64),private,parameter :: t_coords(nn_el) = [-1.0_real64, -1.0_real64, -1.0_real64,&
                                    -1.0_real64,  1.0_real64,  1.0_real64, 1.0_real64,  1.0_real64]
    !> 单元积分点数
    integer,private,parameter :: ngauss=8
    !> 单元高斯积分的等参坐标r数组
    real(real64),private,parameter :: gauss_r_coords(nn_el) = (sqrt(1.0_real64/3.0_real64))*[-1.0_real64,&
    1.0_real64,  1.0_real64, -1.0_real64, -1.0_real64,  1.0_real64, 1.0_real64, -1.0_real64]
    !> 单元高斯积分的等参坐标s数组
    real(real64),private,parameter :: gauss_s_coords(nn_el) = (sqrt(1.0_real64/3.0_real64))*[-1.0_real64,&
    -1.0_real64,  1.0_real64,  1.0_real64, -1.0_real64, -1.0_real64, 1.0_real64,  1.0_real64]
    !> 单元高斯积分的等参坐标t数组
    real(real64),private,parameter :: gauss_t_coords(nn_el) = (sqrt(1.0_real64/3.0_real64))*[-1.0_real64,&
    -1.0_real64, -1.0_real64, -1.0_real64,  1.0_real64,  1.0_real64, 1.0_real64,  1.0_real64]
    !> 单元高斯积分的权重数组
    real(real64),private,parameter :: gauss_weights(nn_el) = [1.0_real64, 1.0_real64, 1.0_real64, 1.0_real64,&
    1.0_real64, 1.0_real64, 1.0_real64, 1.0_real64]

    contains

    !> 计算C3D8单元刚度矩阵(b-bar)
    subroutine c3d8_ke(e, nu, el_ncoords,ke)
        real(real64), intent(in) :: e, nu, el_ncoords(3,nn_el)
        real(real64), intent(out) :: ke(eldofs,eldofs)
        call get_bbar_ke(e, nu, el_ncoords, ke)
    end subroutine c3d8_ke

    !> 计算C3D8单元的集中质量矩阵
    subroutine c3d8_me(rho, el_ncoords, gme)
        real(real64), intent(in) :: rho,el_ncoords(3,nn_el)
        real(real64), intent(out) ::  gme(eldofs,eldofs)
        call get_lumped_me(rho, el_ncoords, gme)
    end subroutine c3d8_me

    !> 计算C3D8单元刚度矩阵(wilson element)
    subroutine c3d8i_ke(e, nu, el_ncoords,ke)
        real(real64), intent(in) :: e, nu, el_ncoords(3,nn_el)
        real(real64), intent(out) :: ke(eldofs,eldofs)
        call get_wilson_ke(e, nu, el_ncoords, ke)
    end subroutine c3d8i_ke

    !> 计算C3D8I单元的集中质量矩阵
    subroutine c3d8i_me(rho, el_ncoords, gme)
        real(real64), intent(in) :: rho,el_ncoords(3,nn_el)
        real(real64), intent(out) ::  gme(eldofs,eldofs)
        call get_lumped_me(rho, el_ncoords, gme)
    end subroutine c3d8i_me
! ------------------------------------------
!   以下为 模块内部使用的子程序和函数
! -------------------------------------------
    !> 计算第i个节点的第dofth自由度的全局自由度
    function get_dof_ind(i,dofth) result(dof_ind)
        integer,intent(in):: i,dofth
        integer::dof_ind
        select case(dofth)
            case(1)
                dof_ind = i*3-2
            case(2)
                dof_ind = i*3-1
            case(3)
                dof_ind = i*3
        end select
    end function get_dof_ind

    !> 计算节点的自由度索引数组(3)
    function get_ndof_inds(i) result(retval)
        integer, intent(in) :: i
        integer :: retval(3)
        retval(1) = nndofs*i-2
        retval(2) = nndofs*i-1
        retval(3) = nndofs*i
    end function get_ndof_inds

    !> 计算单元的全局自由度索引数组
    subroutine get_eldof_inds(el_conn,el_inds)
        implicit none
        integer, intent(in) :: el_conn(nn_el)
        integer, intent(out) ::  el_inds(eldofs)
        integer:: i,j
        el_inds(:)=0
        do i=1,nn_el
            j=el_conn(i)
            el_inds(3*i-2)  = 3*j-2
            el_inds(3*i-1)  = 3*j-1
            el_inds(3*i)    = 3*j
        end do
    end subroutine get_eldof_inds
    
    !> 计算各向同性弹性矩阵 D (sigma = D * epsilon)，单精度
    function get_d_matrix(E, nu) result(d)
        real(real64), intent(in)  :: E           ! 弹性模量
        real(real64), intent(in)  :: nu          ! 泊松比
        real(real64) :: d(6,6)  ! 输出: 6x6 弹性矩阵
        ! 初始化矩阵为零
        d = 0.0_real64; 
        ! 对角项 (1-nu)
        d(1,1) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * (1.0_real64 - nu)
        d(2,2) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * (1.0_real64 - nu)
        d(3,3) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * (1.0_real64 - nu)
        ! 耦合项 (nu)
        d(1,2) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * nu
        d(1,3) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * nu
        d(2,3) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * nu
        d(2,1) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * nu
        d(3,1) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * nu
        d(3,2) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * nu
        ! 剪切项 (1-2*nu)*0.5
        d(4,4) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * (1.0_real64 - 2.0_real64 * nu) * 0.5_real64
        d(5,5) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * (1.0_real64 - 2.0_real64 * nu) * 0.5_real64
        d(6,6) = (E / ((1.0_real64 + nu) * (1.0_real64 - 2.0_real64 * nu))) * (1.0_real64 - 2.0_real64 * nu) * 0.5_real64
    end function get_d_matrix

    !> 返回形函数Ni在(r,s,t)的值
    function get_ni_val(i,r,s,t) result(retval)
        integer, intent(in) :: i
        real(real64), intent(in) :: r,s,t
        real(real64) :: retval
        real(real64) :: r_i,s_i,t_i
        r_i = r_coords(i)
        s_i = s_coords(i)
        t_i = t_coords(i)
        retval = 0.125_real64*(1.0_real64+r_i*r)*(1.0_real64+s_i*s)*(1.0_real64+t_i*t)
    end function get_ni_val
    
    !> 返回i节点形函数Ni对r的偏导
    function get_dni_dr(i,s,t) result(retval)
        integer, intent(in) :: i
        real(real64), intent(in) :: s,t
        real(real64) :: retval
        real(real64) :: r_i,s_i,t_i
        r_i = r_coords(i)
        s_i = s_coords(i)
        t_i = t_coords(i)
        retval = 0.125_real64*r_i*(1.0_real64+s_i*s)*(1.0_real64+t_i*t)
    end function get_dni_dr

    !> 返回i节点形函数Ni对s的偏导
    function get_dni_ds(i,r,t) result(retval)
        integer, intent(in) :: i
        real(real64), intent(in) :: r,t
        real(real64) :: retval
        real(real64) :: r_i,s_i,t_i
        r_i = r_coords(i)
        s_i = s_coords(i)
        t_i = t_coords(i)
        retval = 0.125_real64*s_i*(1.0_real64+r_i*r)*(1.0_real64+t_i*t)
    end function get_dni_ds

    !> 返回i节点形函数Ni对t的偏导
    function get_dni_dt(i,r,s) result(retval)
        integer, intent(in) :: i
        real(real64), intent(in) :: r,s
        real(real64) :: retval
        real(real64) :: r_i,s_i,t_i
        r_i = r_coords(i)
        s_i = s_coords(i)
        t_i = t_coords(i)
        retval = 0.125_real64*t_i*(1.0_real64+r_i*r)*(1.0_real64+s_i*s)
    end function get_dni_dt
    
    !> 返回单元在(r,s,t)的形函数矩阵
    function get_nmat(r,s,t) result(nmat)
        real(real64), intent(in) :: r,s,t
        real(real64) :: nmat(nndofs,eldofs)
        integer :: i
        nmat(:,:)=0.0_real64
        do i=1,nn_el
            nmat(1,3*i-2) = get_ni_val(i,r,s,t)
            nmat(2,3*i-1) = get_ni_val(i,r,s,t)
            nmat(3,3*i)   = get_ni_val(i,r,s,t)
        end do
    end function get_nmat

    !> 计算(r,s,t)处的Jacbian矩阵(3x3)
    function get_jacbmat(r,s,t,el_ncoords) result(jac)
        real(real64), intent(in) :: r,s,t,el_ncoords(3,nn_el)
        real(real64) :: jac(3,3)
        integer :: i
        real(real64) :: ncoords(nn_el,3), dn_p(3,8) ! 节点坐标矩阵(8,3)和形函数导数矩阵(3,8)
        ncoords=transpose(el_ncoords)
        dn_p(:,:)=0.0_real64
        do i=1,nn_el
            dn_p(1,i)=get_dni_dr(i,s,t)
            dn_p(2,i)=get_dni_ds(i,r,t)
            dn_p(3,i)=get_dni_dt(i,r,s)
        end do
        jac = matmul(dn_p,ncoords)    ! 计算Jacobian矩阵
    end function get_jacbmat

    !> 计算(r,s,t)处的Jacbian矩阵的行列式
    function get_jacb_det(r,s,t,el_ncoords) result(det)
        use matrix_utils_mod, only: det3rank
        implicit none
        real(real64), intent(in) :: r,s,t
        real(real64), intent(in) :: el_ncoords(3,nn_el)
        real(real64) :: det,jac(3,3)
        integer :: info

        jac = get_jacbmat(r,s,t,el_ncoords)
        call det3rank(jac,det,info)
        if (info /= 0) then
            stop 'Error in function get_jacb_det'
        end if
    end function get_jacb_det

! -------------------------------------------------------------
!   use b-bar modified technique to solve volume locking
! -------------------------------------------------------------
    !> b-bar修正应变矩阵b_dev, b_dil
    subroutine get_bbar_bmat(r,s,t,el_ncoords,b_dev,b_dil)
        !> el_ncoords : 单元节点坐标数组(3,8)
        use matrix_utils_mod
        implicit none
        real(real64), intent(in)  :: r,s,t,el_ncoords(3,nn_el)
        real(real64), intent(out) :: b_dev(6,eldofs), b_dil(6,eldofs)
        real(real64):: t_mat(6,6),jac(3,3),jac_inv(3,3),dn_dris(3),dn_dxis(3)&
        ,bi_standard(6,3),bi_dil(6,3),bi_dev(6,3)
        integer:: i,info
        
        b_dev(:,:) =0.0_real64
        b_dil(:,:) =0.0_real64
        ! 赋值t矩阵
        t_mat(:,:) = 0.0_real64
        do i = 1, 3
            t_mat(i,i) = 1.0_real64
        end do
        t_mat(4,5) = 1.0_real64
        t_mat(5,6) = 1.0_real64
        t_mat(6,4) = 1.0_real64
        
        ! get jacabian mat at (r,s,t)
        jac = get_jacbmat(r,s,t,el_ncoords)
        call invert_matrix_d(jac, jac_inv, info) ! 计算逆矩阵
        if (info /= 0) then
            print *, 'Error in subroutine get_bbar_bmat'
            stop
        end if

        ! 计算b_dev和b_dil
        do i = 1, 8
            dn_dris(1) = get_dni_dr(i,s,t)
            dn_dris(2) = get_dni_ds(i,r,t)
            dn_dris(3) = get_dni_dt(i,r,s)
            ! 形函数对x,y,z的偏导数组: Nx,Ny,Nz
            dn_dxis = matmul(jac_inv,dn_dris)   
            ! 用Nx,Ny,Nz给Bi赋值
            bi_standard(:,:) = 0.0_real64
            bi_standard(1,1) = dn_dxis(1)
            bi_standard(2,2) = dn_dxis(2)
            bi_standard(3,3) = dn_dxis(3)
            bi_standard(4,1) = dn_dxis(2)
            bi_standard(4,2) = dn_dxis(1)
            bi_standard(5,2) = dn_dxis(3)
            bi_standard(5,3) = dn_dxis(2)
            bi_standard(6,1) = dn_dxis(3)
            bi_standard(6,3) = dn_dxis(1)
            bi_standard = matmul(t_mat,bi_standard)

            bi_dil(:,:) = 0.0_real64
            bi_dil(1,1) =dn_dxis(1)/3.0_real64
            bi_dil(2,1) =dn_dxis(1)/3.0_real64
            bi_dil(3,1) =dn_dxis(1)/3.0_real64
            bi_dil(1,2) =dn_dxis(2)/3.0_real64
            bi_dil(2,2) =dn_dxis(2)/3.0_real64
            bi_dil(3,2) =dn_dxis(2)/3.0_real64
            bi_dil(1,3) =dn_dxis(3)/3.0_real64
            bi_dil(2,3) =dn_dxis(3)/3.0_real64
            bi_dil(3,3) =dn_dxis(3)/3.0_real64

            bi_dev = bi_standard-bi_dil
            b_dil(:,3*i-2:3*i) = bi_dil
            b_dev(:,3*i-2:3*i) = bi_dev
        end do
    end subroutine get_bbar_bmat

    subroutine get_bbar_ke(e,nu,el_ncoords,ke)
        real(real64), intent(in)  :: el_ncoords(3,nn_el),e,nu
        real(real64), intent(out) :: ke(eldofs,eldofs)
        real(real64)              :: ri,si,ti,wi,detjac
        real(real64)              :: dmat(6,6),b_dev(6,eldofs),b_dil(6,eldofs)
        integer                   :: i
        
        ke(:,:)    = 0.0_real64
        b_dev(:,:) = 0.0_real64
        b_dil(:,:) = 0.0_real64
        dmat = get_d_matrix(e,nu)
        ! 对b_dev进行全积分
        do i=1,ngauss
            ri = gauss_r_coords(i)
            si = gauss_s_coords(i)
            ti = gauss_t_coords(i)
            wi = gauss_weights(i)
            detjac = get_jacb_det(ri,si,ti,el_ncoords)
            call get_bbar_bmat(ri,si,ti,el_ncoords,b_dev,b_dil)
            ke=ke+wi*detjac*matmul(transpose(b_dev),matmul(dmat,b_dev))
        end do

        ! 对b_dil进行缩减积分
        ri = 0.0_real64
        si = 0.0_real64
        ti = 0.0_real64 
        wi = 8.0_real64
        call get_bbar_bmat(ri,si,ti,el_ncoords,b_dev,b_dil)
        detjac = get_jacb_det(ri,si,ti,el_ncoords)
        ke=ke+wi*detjac*matmul(transpose(b_dil),matmul(dmat,b_dil))
    end subroutine get_bbar_ke
    
    !> 计算单元的集中质量矩阵: 将单元总质量平均分到8个节点
    subroutine get_lumped_me(rho, el_ncoords, gme)
        real(real64), intent(in)    :: rho,el_ncoords(3,nn_el)
        real(real64), intent(out)   :: gme(eldofs,eldofs)
        real(real64)                :: el_mass,wi,ri,si,ti,detjac
        integer                     :: i
        el_mass  = 0.0_real64
        gme(:,:) = 0.0_real64
        do i=1,ngauss
            ri      = gauss_r_coords(i)
            si      = gauss_s_coords(i)
            ti      = gauss_t_coords(i)
            wi      = gauss_weights(i)
            detjac  = get_jacb_det(ri,si,ti,el_ncoords)
            el_mass = el_mass+wi*detjac*rho
        end do
        do i = 1,eldofs
            gme(i,i)=0.125_real64*el_mass
        end do
    end subroutine get_lumped_me
! -------------------------------------------------------------
!   use wilson imcompatible displacement mode to solve shear locking
!   from : qinguoming.xyz
! -------------------------------------------------------------
    !> 计算节点i的附加自由度形函数在(r,s,t)处的值
    function get_mi_val(i,r,s,t) result(retval)
        real(real64), intent(in) :: r,s,t
        real(real64) :: retval
        integer, intent(in) :: i
        select case (i)
            case (1)
                retval = 1.0_real64-r*r
            case (2)
                retval = 1.0_real64-s*s
            case (3)
                retval = 1.0_real64-t*t
        end select 
    end function get_mi_val

    !> 返回i节点形函数Mi对r的偏导在(r,s,t)处的值
    function get_dmi_dr(i,r) result(retval)
        real(real64), intent(in) :: r
        integer, intent(in)::i
        real(real64) :: retval
        select case (i)
            case(1)
                retval = -2.0_real64*r
            case(2)
                retval = 0.0_real64
            case(3)
                retval = 0.0_real64
        end select
    end function get_dmi_dr
    
    !> 返回i节点形函数Mi对s的偏导在(r,s,t)处的值
    function get_dmi_ds(i,s) result(retval)
        real(real64), intent(in) :: s
        integer, intent(in)::i
        real(real64) :: retval
        select case (i)
            case(1)
                retval = 0.0_real64
            case(2)
                retval = -2.0_real64*s
            case(3)
                retval = 0.0_real64
        end select
    end function get_dmi_ds

    !> 返回i节点形函数Mi对t的偏导在(r,s,t)处的值
    function get_dmi_dt(i,t) result(retval)
        real(real64), intent(in) :: t
        integer, intent(in)::i
        real(real64) :: retval
        select case (i)
            case(1)
                retval = 0.0_real64
            case(2)
                retval = 0.0_real64
            case(3)
                retval = -2.0_real64*t
        end select
    end function get_dmi_dt

    !> 计算(r,s,t)处的附加自由度的应变矩阵bimat
    function get_wilson_bmat(r,s,t,el_ncoords) result(wb)
        use matrix_utils_mod
        implicit none
        real(real64), intent(in)    :: r,s,t,el_ncoords(3,nn_el)
        real(real64)                :: wb(6,9)
        real(real64)                :: jac(3,3),jac_inv(3,3),dm_dris(3),dm_dxis(3)
        integer                     :: i,info

        !通过jacbian矩阵计算Mi对x,y,z的偏导数组
        jac = get_jacbmat(r,s,t,el_ncoords)
        call invert_matrix_d(jac, jac_inv, info) ! 计算逆矩阵
        if (info /= 0) then
            print *, 'Error in subroutine get_wilson_bmat'
        end if
        ! 组装wilson B mat
        wb(:,:) = 0.0_real64
        do i=1,3
            !计算Mi对r,s,t的偏导数组
            dm_dris(1) = get_dmi_dr(i,r)
            dm_dris(2) = get_dmi_ds(i,s)
            dm_dris(3) = get_dmi_dt(i,t)
            ! 计算Mi对x,y,z的偏导数组
            dm_dxis = matmul(jac_inv,dm_dris)
            ! 用Mx,My,Mz给B赋值
            wb(1,3*i-2) = dm_dxis(1)
            wb(2,3*i-1) = dm_dxis(2)
            wb(3,3*i)   = dm_dxis(3)
            wb(4,3*i-2) = dm_dxis(2)
            wb(4,3*i-1) = dm_dxis(1)
            wb(5,3*i-1) = dm_dxis(3)
            wb(5,3*i)   = dm_dxis(2)
            wb(6,3*i-2) = dm_dxis(3)
            wb(6,3*i)   = dm_dxis(1)
        end do
    end function get_wilson_bmat

    !>计算(r,s,t)处, C3D8原始单元的应变矩阵B
    function get_hex_bmat(r,s,t,el_ncoords) result(bmat)
        use matrix_utils_mod
        implicit none
        real(real64),intent(in) :: r,s,t,el_ncoords(3,nn_el)
        real(real64) :: bmat(6,eldofs)
        real(real64)::jac(3,3),jac_inv(3,3),dn_dris(3),dn_dxis(3)
        integer::i,info
        
        jac = get_jacbmat(r,s,t,el_ncoords)
        call invert_matrix_d(jac, jac_inv, info) ! 计算逆矩阵
        if (info /= 0) then
            print *, 'error in function : get_hex_bmat'
            stop
        end if
        do i=1,nn_el
            dn_dris(1) = get_dni_dr(i,s,t)
            dn_dris(2) = get_dni_ds(i,r,t)
            dn_dris(3) = get_dni_dt(i,r,s)
            dn_dxis = matmul(jac_inv,dn_dris)
            bmat(1,3*i-2) = dn_dxis(1)
            bmat(2,3*i-1) = dn_dxis(2)
            bmat(3,3*i)   = dn_dxis(3)
            bmat(4,3*i-2) = dn_dxis(2)
            bmat(4,3*i-1) = dn_dxis(1)
            bmat(5,3*i-1) = dn_dxis(3)
            bmat(5,3*i)   = dn_dxis(2)
            bmat(6,3*i-2) = dn_dxis(3)
            bmat(6,3*i) = dn_dxis(1)
        end do
    end function get_hex_bmat

    !> 计算根据wilson非协调单元修正的单元刚度矩阵
    subroutine get_wilson_ke(e,nu,el_ncoords,ke)
        use matrix_utils_mod
        implicit none
        real(real64), intent(in) :: el_ncoords(3,nn_el),e,nu
        real(real64), intent(out) :: ke(eldofs,eldofs)
        integer::i,info
        real(real64)::ri,si,ti,wi,detjac
        real(real64)::dmat(6,6),kcc(eldofs,eldofs),kii(9,9),kii_inv(9,9),&
                kci(eldofs,9),kic(9,eldofs),bcmat(6,eldofs),bimat(6,9)
        ! 初始化
        kcc(:,:)    = 0.0_real64
        kii(:,:)    = 0.0_real64
        kci(:,:)    = 0.0_real64
        kic(:,:)    = 0.0_real64
        bcmat(:,:)  = 0.0_real64
        bimat(:,:)  = 0.0_real64
        dmat = get_d_matrix(e,nu)
        ! 先计算kcc,kii,kci,kic矩阵
        do i=1,ngauss
            ri = gauss_r_coords(i)
            si = gauss_s_coords(i)
            ti = gauss_t_coords(i)
            wi = gauss_weights(i)
            detjac = get_jacb_det(ri,si,ti,el_ncoords)
            ! 计算b矩阵
            bcmat = get_hex_bmat(ri,si,ti,el_ncoords)
            bimat = get_wilson_bmat(ri,si,ti,el_ncoords)
            kcc = kcc + wi*detjac*matmul(transpose(bcmat),matmul(dmat,bcmat))
            kii = kii + wi*detjac*matmul(transpose(bimat),matmul(dmat,bimat))
            kci = kci + wi*detjac*matmul(transpose(bcmat),matmul(dmat,bimat))
            kic = kic + wi*detjac*matmul(transpose(bimat),matmul(dmat,bcmat))
        end do

        ! 静力凝聚
        call invert_matrix_d(kii,kii_inv,info)
        if ( info /= 0 ) then
            print *, 'Error in subroutine get_wilson_ke'
            stop
        end if
        ke = kcc - matmul(kci,matmul(kii_inv,kic))
    end subroutine get_wilson_ke
end module solid_mod